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Upon sensory stimulation, primary cortical areas readily engage in narrow-band rhythmic 
activity between 30 and 90 Hz, the so-called gamma oscillations. Here we show that, 
when embedded in a balanced network, type-l excitable neurons entrained to the 
collective rhythm show a discontinuity in their firing-rates between a slow and a fast 
spiking mode. This jump in the spiking frequencies is characteristic to type II neurons, 
but is not present in the frequency-current curve [f-l curve) of isolated type I neurons. 
Therefore, this rate bimodality arises as an emerging network property in type I population 
models. We have studied the mechanisms underlying the generation of these two firing 
modes, in order to reproduce the spiking activity of in vivo cortical recordings, which is 
known to be highly irregular and sparse. We have also analyzed the relation between 
afferent inputs and the single unit activity, and between the latter and the local field 
potential (LFP) phase, in order to establish how the collective dynamics modulates the 
spiking activity of the individual neurons. Our results reveal that the inhibitory-excitatory 
balance allows two encoding mechanisms, for input rate variations and LFP phase, to 
coexist within the network. 
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1. INTRODUCTION 

The type and distribution of ionic channels across the membrane 
of a neuron determine its excitability behavior. This behavior 
can be experimentally tested by injecting a pulsed, constant, or 
ramp current. In the latter case one observes how the membrane 
potential changes as the input amplitude increases, eventually 
switching from a resting to an oscillating (tonic firing) state. 
Mathematically, this transition is described as a bifurcation of the 
membrane potential from a stable fixed point to a periodic orbit. 

Depending on how the firing frequency within the oscil- 
latory regime behaves near the bifurcation, neural excitability 
can be classified into type I or type II (Rinzel and Ermentrout, 
1989). In type I excitability, the frequency increases continu- 
ously from zero as the bifurcation is crossed. Type II neurons, 
on the other hand, exhibit a discontinuous jump in frequency 
as the tonic regime is entered, and their range of firing fre- 
quencies is quite narrow compared to type I neurons, which can 
achieve arbitrarily low frequencies (Izhikevich, 2007). These dif- 
ferent behaviors are associated with distinct bifurcations. In type 
I excitability, the periodic orbit emerges due to the collision of a 
stable fixed point (the resting state) and an unstable equilibrium 
point, occurring on top of an invariant circle (SNIC bifurca- 
tion). Type II excitability, on the other hand, can arise in three 
different ways: via a subcritical or a supercritical Hopf bifurca- 
tion, or through a saddle-node bifurcation outside the invariant 
circle (Rue and Garcia-Ojalvo, 2011). The integrate-and-fire and 



conductance-based models used in the literature to describe cor- 
tical networks are usually of type I (Wang and Buzsaki, 1996; 
El Boustani etal., 2007). 

The phase response can also be used as a criterion to dis- 
tinguish between excitability classes: when operating in a tonic 
regime, type I neurons always advance their phase (defined with 
respect to their spiking period) when perturbed by a brief depo- 
larizing pulse, while type II neurons can either advance or delay 
their phase depending on the instant of perturbation relative to 
the period of oscillation (Hansel et al., 1995). Within a network, 
the input pulses come from other neurons, and therefore synchro- 
nized periodic spiking might lead to lower or higher frequency 
rhythms depending on the type of neuronal excitability. However, 
high spike synchrony and periodic firing patterns are rarely seen 
in vivo; in fact the coefficient of variation of interspike intervals 
in the cortical neurons is typically larger than 1 (Shinomoto et al, 
1999, 2005). Consequently, the rhythm of the neural population 
must emerge from the recurrent synaptic activity in the network, 
rather than from spike-to-spike synchrony. It has been experi- 
mentally observed (Softky and Koch, 1993) that the high synaptic 
bombardment acting upon cortical neurons is far from being 
constant, leading to stochastic fluctuations that affect these neu- 
rons in vivo. Thus, neither the single-neuron frequency- current 
curve if -I curve) nor the phase response can perfectly describe 
the behavior of single neurons when they are embedded in a 
network. 
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Cortical oscillations thus arise as a collective phenomenon, 
which does not require individual-neuron firings to be oscillatory 
themselves. These global oscillations are discernible in averaged 
population activities such as the local field potential (LFP), whose 
troughs (or peaks, depending on the relative position of the 
recording electrode and the generating current dipoles) corre- 
spond to the minima of the synaptic inhibitory flow. In these 
temporal windows neurons are more likely to spike, produc- 
ing an increase in the excitatory synaptic current, followed by a 
new burst of inhibition. These oscillations, at frequencies in the 
beta (12-30 Hz) and gamma (30-90 Hz) ranges, are experimen- 
tally seen in the cortex upon sensory stimulation (Buzsaki and 
Draguhn, 2004). Specifically, gamma-band synchronization has 
received special consideration, as it is hypothesized to be a mecha- 
nism for the dynamic generation of functional cell assemblies and 
for the flexible communication between brain areas [for a review, 
see Singer (1999) and Fries (2009)]. 

Here we study the generation of gamma-band (30-90 Hz) 
oscillations that occur without spike -to -spike synchrony 
[known as synchronous irregular dynamics (van Vreeswijk and 
Sompolinsky, 1996; Brunei, 2000; Hansel and Mato, 2003)], and 
how these rhythms affect single unit activity. Specifically, we 
consider networks of type I excitatory and inhibitory neurons 
described by conductance-based models and stimulated by 
slowly varying excitatory inputs, and study the emergent firing 
patterns of the neurons once embedded in a balanced network. 
Our results show that the synchronous irregular state giving rise 
to gamma-band oscillations is composed of two firing modes: a 
high-frequency (or bursty) regime and a low-frequency regime 
separated by a gap of quasi-forbidden firing frequencies. We also 
study the relationship between this firing behavior and the input: 
the fast spiking mode encodes for the stimulus rate, whereas 
the existence of a slow spiking mode allows for a phase code to 
appear. 

2. METHODS 

2.1. DESCRIPTION OF THE NEURONAL MODELS 

The dynamical equation for the neuronal membrane voltage is 
given by a conductance-based model: 



cj^- = -g K n\V-V K ) -gN a m 3 h(V-V Na ) 
at 

-gdv- v L ) + i syn 



(i) 



where gK , gNa> and gi are the maximum conductances for the 
potassium, sodium and the leak current, respectively, and J syn 
is the synaptic current coming from the neighboring neurons 
impinging on one neuronal cell. The dynamics of the potas- 
sium and sodium channels is represented by the time-varying 
probabilities that a channel is open: 



dx 
It 



[a x (V)(l-x)-V x (V)x], 



where x stands for n in the case of the potassium current, and 
for m and h in the case of the sodium current. a(V) and f5(V) 
are voltage-dependent rate constants, and (f> is the temperature 



factor, defined by c|> = 3( T ~ 6 - 3 )/ 10 ; where T is measured in degrees 
Celsius. 

The parameter values used throughout this study are those 
of Ref. Gutfreund et al. (1995): g K = 4.74 |xS, g Na = 12.5 |xS, 
and gi = 0.025 |xS. The reversal potentials of the different chan- 
nels are V K = -80 mV, V Na = 40 mV, and V L = -65 mV, and 
the membrane capacitance is C m = 0.25 nF (0.125 nF) for the 
excitatory (inhibitory)neurons. The leak conductance defines 
an effective membrane time constant for the isolated neuron 
according to the expression x = C m /gi, which is taken to be 
10 ms for the excitatory neurons and 5 ms for the inhibitory 
neurons. The temperature factor (j> is set to 21, which corre- 
sponds to T = 34°C. These parameter values lead to type I 
excitability. 

The rate functions ot and f5 for each gating variable are: 



a„(V) = 0.01- 



V + 20 



1 _ e -(V+20)/10 

P„(V) = 0.125 e - (V+30) / 80 



for the gating variable n, 



a m (V) = 0.1- 



V+ 16 



1 _ e -(V+l6)/W 

(J m (V) = 4 e -< v + 41 >/i8 



for the gating variable m, and 



O fc (V) = 0.07 6 - (V '+ 30) / 20 
1.0 



l + g-V/10 



for the gating variable h. Due to the rapid activation of m we 
replace it by its steady-state value moo = a +g ■ 

2.2. DESCRIPTION OF THE NETWORK MODEL 

We consider a network composed of 2000 neurons, 80% of 
which are excitatory while the remaining 20% are inhibitory 
(Soriano et al, 2008). All connections between cells are chemi- 
cal synapses — no gap junctions are considered — and each neuron 
connects randomly with 200 other neurons, on average, belong- 
ing to both populations. Therefore, no architecture is imposed 
on the connectivity. We have also introduced a synaptic delay 
in the transmission of the action potential between neighbor- 
ing neurons, taken from a gamma distribution of mean 2 ms and 
variance 4 ms 2 . The synaptic current is described using again a 
conductance-based formalism: 



hyn = gsyn(t)(V{t) - E syn ), 

where g syn is the synaptic conductance and £ syn is the reversal 
potential of the synapse. For £ syn greater than the resting poten- 
tial Vj-est the synapse is depolarizing, i.e., excitatory, otherwise it is 
hyperpolarizing, i.e., inhibitory. We consider two temporal time 
constants, Xd and x, (decay and rise synaptic time, respectively, 
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Table 1 | Synaptic time constants and reversal synaptic potential 
values. 



Synapse 



x r (ms) 



x d (ms) 



(mV) 



AM PA 
GABA 



0.5 
2 



0 

-70 



The equilibrium potential at l syn = 0 for Equation (11 is V rest = —65mV. 



Table 2 | Synaptic conductances, g' svn , for all the possible 
connections. 



Synapse 

GABA 

Recurrent AM PA 
External AM PA 



Conductance on Conductance on 

inhibitory neurons (nS) excitatory neurons (nS) 



240 
2.5 
3.2 



240 
2.5 
3.2 



see Table 1), for the dynamics of the synaptic conductance, which 
is calculated by 



gsyn(f) 



bsyn 



where g ' , shown in Table 2, is tuned in order to obtain a balance 
between excitation and inhibition, given the/-/ relation. 

We have chosen the maximal conductances, g' SYn , to maintain 
the postsynaptic potential (PSP) amplitudes within physiological 
ranges: the excitatory PSP in the range from 0.42 to 0.83 mV, and 
the inhibitory PSP from 1.54 to 1.88 mV. 

Additionally, all neurons receive an heterogeneous Poisson 
train of excitatory presynaptic potentials with a mean event 
rate that varies following an Ornstein-Uhlenbeck process. This 
incoming external current mimics the direct input from any other 
area external to the network considered here. The instantaneous 
rate, ~k(t), of the external excitatory train of spikes is generated 
according to an Ornstein-Uhlenbeck process as considered in 
Mazzoni et al. (2008): 



dX 
di 



-X(f)+a(0 



!](*)> 



where o(t) is the standard deviation of the noisy process and is 
set to 0.6spikes/s. x is set to 16 ms, leading to a 1/f power spec- 
trum for the X time series that is flat up to a cut-off frequency 
/ = = 9.9 Hz. r\(t) is a Gaussian white noise. 

Excitatory synapses outnumber inhibitory ones, and yet the 
brain avoids epileptic states because inhibition is able to balance 
excitation, and thus neurons remain below threshold, firing only 
occasionally. Inhibitory neurons have higher firing-rates than 
excitatory neurons for any given input current, as can be seen in 
the/-!" curve of Figure 1. Additionally, GABAergic (inhibitory) 
synapses are stronger than glutamatergic AMPA (excitatory) 
synapses to compensate their relative small number (Markram 
etal.,2004). 



2.3. NUMERICAL SIMULATIONS 

The model has been integrated using the Heun algorithm (Garcia- 
Ojalvo and Sancho, 1999), with a time step of 0.05 ms. All the 
simulations represent 3 s of activity and the connectivity, initial 
conditions, and noise realization were varied from trial to trial. 

2.4. MODEL OF LFP 

We quantify the activity of the network in different ways. At the 
single-neuron level we consider the instantaneous firing-rate as 
a measure of the individual spiking dynamics. At the population 
level we use two observables: the time-resolved average firing-rate 
of the whole neuronal population (defined as the total number of 
spikes per unit time in the population divided by the number of 
neurons), and the LFP, computed as the sum of the absolute val- 
ues of the excitatory and inhibitory synaptic currents acting upon 
the excitatory neurons, averaged over this population (Mazzoni 
et al, 2008; Buzsaki et al, 2012): 



LFP = R e <|/ampaI + IXgabaI), 



(2) 



Here (• • • ) denotes an average over all excitatory neurons (Berens 
et al., 2010). The term -Tampa accounts for both the external 
excitatory heterogeneous Poisson spike train and the recurrent 
excitatory synaptic current due to network connectivity, while 
Jgaba corresponds to the recurrent inhibitory synaptic current. 
R e represents the resistance of a typical electrode used for extracel- 
lular measurements, here chosen to be 1 MQ,. In Figure 9 the LFP 
was filtered with a 4th order Butterworth bandpass filter using 
MATLAB function filter.m. 

2.5. SPIKE TRIGGERED AVERAGE 

We have calculated the spike triggered average (STA) of the LFP 
and of the inhibitory synaptic current impinging on the neurons. 
For the considered spikes we have registered these signals during 
a time window starting 50 ms prior to the spike and ending 20 ms 
after it, and computed the mean across the total number of action 
potentials. 

2.6. COMPUTATION OF POWER SPECTRA 

The LFP power spectrum was estimated using the multitaper 
method (Thomson, 1982) commonly used to reduce the variance 
of the spectra of recorded signals, which are usually very noisy. 
This estimator was implemented in Chronux 2.10 (Bokil et al., 
2010). The multitapered power spectrum, S(/), is the average of 
the power spectrum of the LFP signal multiplied by K orthogonal 
Slepian functions (in our case K = 5), and further averaged over 
N = 20 trials: 



S(f) 



N N ( K \ 

n=l n = l \ k=l I 



(3) 



Here LFP nJ t(/) is the discrete Fourier transform of the LFP(f) sig- 
nal of the «-th trial, multiplied by the /c-th Slepian function (or 
taper). We have considered data segments within a 500-ms sliding 
time window with an overlap of 50 ms, padded with zeros up to a 
length of 512 in order to obtain an increased sampling rate in the 
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frequency domain. The resolution bandwidth is thus ±6 Hz. The 
firing-rate power spectra are also obtained by the multitaper algo- 
rithm, with the same sliding time window, overlap, and padding. 
In this case the average instantaneous firing-rate is obtained from 
an histogram of the spiking times, with a 1-ms bin. All histograms 
and power spectra are averaged over 20 trials. 

3. RESULTS 

An adequate neural network model should reproduce several 
experimentally observed features of in vivo cortical oscillations: a 
prominent peak of the LFP power spectrum in the gamma-band 
resulting from external stimulation, irregular individual firing at 
frequencies lower than the gamma rhythm, and partial phase 
locking of individual spikes to the gamma cycle. We show in 
what follows that the neural network described above is able to 
reproduce these features. 

3.1. FIRING-RATE OF NETWORK-EMBEDDED TYPE I NEURONS 

In this section we study how the network shapes the firing pat- 
terns of the single neuron components. Before studying the firing 
modes of the cells when embedded in a network, we first char- 
acterize their behavior as single cells. When isolated, a neuron 
displays a characteristic transition to the tonic firing regime as the 
input current is gradually increased (solid lines in Figure 1). Type 
I integrator cells can achieve arbitrary low spiking frequencies. 
This behavior is a consequence of the saddle-node bifurcation 
on an invariant circle that this type of neuron undergoes with 
increasing injected current. For the same current, inhibitory cells 
(with a smaller membrane time constant) fire with a higher 
rate. However, both excitatory and inhibitory cells have the same 
spiking threshold. 

The behavior of the firing-rate under constant injection cur- 
rent, shown in Figure 1 A for adiabatically increasing values of the 
current, is not suited for characterizing the behavior of a neu- 
ron within a network. When embedded in a population, neurons 



A 




0.2 0.4 0.6 0.8 1 1.2 1.4 

input current (nA) 

FIGURE 1 | (A) f-l curve of type I conductance-based neurons in isolation, 
obtained with the XPPAUT software (Ermentrout, 2002). (B) Equivalent f-l 
curve for neurons embedded in the network (circles) representing the 
instantaneous firing-rate versus the net synaptic current averaged over the 



are subject to synaptic input that fluctuates in time. Even when 
the firing of the presynaptic neurons projecting to a given postsy- 
naptic cell is uncorrected between them, the total input current 
is Gaussian distributed (Roxin et al., 2011). This input is the 
sum of the synaptic currents coming from all possible sources, 
which include recurrent excitation and inhibition (generated by 
the network itself), and afferent pathways terminating onto the 
population (represented in the model by a Poisson-distributed 
external train of spikes). The net synaptic current impinging 
on a neuron has a nearly stationary average value and fluctu- 
ates rhythmically around it (due to the network-induced gamma 
oscillations), crossing randomly the threshold and leading to an 
irregular individual firing-rate. Here, we capitalize on the nat- 
ural variability of such recurrent and external synaptic drive to 
explore the effective f-I curve of neurons once embedded in a 
balanced recurrent network. Figure IB (circles) represents the 
f-I response curve of type I neurons when embedded within 
a population. Given the above-mentioned fluctuating character 
of the synaptic current received by the neurons within the net- 
work, the f-I response is quantified here in terms of the interval 
between two consecutive spikes (interspike interval, or ISI) as a 
function of the mean synaptic current received by the neuron 
within that interval. This quantity is plotted in Figure IB (cir- 
cles), comparing its behavior with the corresponding/ -I curves of 
the isolated neurons (solid lines). Note, however, the limitations 
of this comparison, as the mean synaptic current over an ISI dif- 
fers from the actual fluctuating input received by a cell, specially 
for long periods. It is precisely in the range of low frequencies 
where both curves differ most. For instance, the firing-rate in the 
resting state region departs from zero close to the bifurcation, 
due to the appearance of noise-induced spikes throughout the 
network. 

Differences between the firing-rate of a network-embedded 
neuron and its characteristic/-/ curve in isolation also exist in the 
tonic regime, as can be seen in Figure IB. We can understand this 



B 




input current (nA) 

corresponding inter-spike interval, for an excitatory (black) and inhibitory (red) 
neuron. In that panel the single-neuron f-l curves are also shown in solid 
lines for comparison. The network is excited by a train of spikes of rate 
8500 spikes/s. 
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from the fact that, due to the synaptic bombardment, neurons 
within a network are persistently perturbed, and consequently 
their firing-rates are altered. The characteristic phase response 
curve (PRC) of type I neurons (Figure 2) reveals a phase advance 
for any perturbation time within a period. Thus, when a type I 
neuron has just spiked and is excited again by a presynaptic poten- 
tial its period is reduced (and thus its firing- rate is increased). 
This explains why the firing-rate of type I neurons embedded in 
a network, and thus continually perturbed by presynaptic acti- 
vation, is above the / -I curve of an isolated neuron with tonic 
firing. 

It is worth noting that the firing-rate discontinuity shown in 
Figure IB (circles) is an emergent network property. Within a 
network, type I neurons display two clearly defined groups of 
inter-spike intervals, in a sort of population-level bimodality that 
does not exist for individual neurons in isolation. 

3.2. NETWORK OSCILLATIONS 

When the induced spiking activity is large enough to percolate 
the network, the internal synaptic currents exhibit a fast cyclic 
behavior, alternating epochs of high excitation followed by high 
inhibition, due to the rhythmic synchronization of GABAergic 
cells via their recurrent connections. The net oscillatory synap- 
tic current leads to a rhythmic behavior in averaged population 
measures such as the LFP and the time-resolved firing-rate of 
the population, as shown in Figures 3A,D, respectively. The LFP 
power spectrum is depicted by solid lines in Figure 3B, reveal- 
ing a frequency peak whose precise position is determined by the 
GABAergic synaptic time constants and the synaptic strengths, as 



well as by the characteristics of the input (Henrie and Shapley, 
2005). For a range of biologically plausible values, the frequency 
of such peak is in the gamma range (30-90 Hz). In particular, for 
the parameters given in Tables 1, 2 (see section 2), the frequency 
peak of the LFP spectrum is located around 45 Hz for an external 
stimulation of mean rate 8500 spikes/s. The low-frequency com- 
ponents are caused by the external incoming spike trains, whose 
mean rate is an Ornstein-Uhlenbeck process. When the effective 
coupling between neurons is low, because either the induced spik- 
ing activity or the synaptic conductances are weak, the collective 
rhythm disappears from the LFP power spectrum, as shown by 
the dashed line in Figure 3B. 

Similarly to the LFP, the global firing-rate exhibits a marked 
peak at ~45Hz in its spectrum (Figure 3E). This rhythmicity 
reflects epochs of synchronization between subsets of neurons. 
Despite the collective rhythmic dynamics (Figure 3D), single 
neurons display strongly irregular individual spiking (Figure 3C), 
firing mostly in a sparse and single-spike mode occasionally 
accompanied by some high-frequency tonic firing, which never- 
theless can be compatible with population rhythmic activity at 
gamma frequencies. 

3.3. FIRING-RATE DISTRIBUTIONS 

As described in the last sections there is a striking difference in 
the spiking activity exhibited by the type I cells once embed- 
ded in a network compared to their intrinsic / -I curve. To better 
resolve the firing mode differences and their origin, we first char- 
acterize the distribution of instantaneous firing-rates, and then 
we describe the relation of firing behavior to the input and LFP 
dynamics. 

To characterize the distribution of instantaneous firing-rates 
of individual neurons within the network described in sec- 
tion 2.2, we compute their histogram for the parameter values 
given in that section, which lead to synchronous irregular fir- 
ing (Figures 3C,D) with global oscillations in the gamma range 
(Figures 3B,E). The distribution of instantaneous firing-rates is 
bimodal (Figure 4A), indicating the prevalence of two firing 
modes, with high-firing-rate events (short ISIs) reflecting bursts 
of spiking activity and short-firing-rate events (long ISIs). This 
effect arises from the network dynamics and cannot be predicted 
from the characteristic f-I curves of type I isolated neurons, 
which are continuous and hence do not forbid any particular 
range of firing-rates (Figure IB; note the gap in the circles, which 
occurs in network-embedded neurons). This bimodal response 
was also observed in a network of LIF neurons (Roxin et al., 201 1), 
when the inhibitory population firing-rate exceeded the excita- 
tory population rate. We will address the origin of this bimodal 
behavior at the end of this section. 

The high-firing-rate peak in Figure 4A arises from the spikes 
of the network-induced bursts (see upper cloud in Figure IB). 
The low-firing-rate events are more frequent than the high- 
firing-rate ones (64% of the interspike intervals of the excitatory 
population across 20 trials lie below ~58.31 spikes/s — see arrows 
in Figure 4A). In other words, more neurons are simultaneously 
found in a silent state or low-firing mode. 

In order to understand the genesis of this bimodal firing 
behavior, we now examine the time evolution of the synaptic 




0 0.2 0.4 0.6 0.8 1 



x/T 

o 

FIGURE 2 | Phase response curve for a type I excitatory and inhibitory 
neuron representing the phase change in the spiking response of the 
neuron to a depolarizing pulse, as a function of the phase of the input 
(normalized to the spiking period T a of the neuron). The excitatory 
(inhibitory) neuron is spiking tonically at an unperturbed period of 8.09 ms 
(6.00 ms), and receives the injection of a pulse with amplitude 1.0 nA and 
duration 0.2 ms. The PRCs, obtained with the XPPAUT software 
(Ermentrout, 2002), are defined as 1 — T(x)/T 0 , where 7~(t) is the period 
after perturbing at time x. 
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FIGURE 3 | Collective oscillations in a type-l-neuron network. (A) LFP 

time trace in a 250 ms interval for an external mean rate of 8500 spikes/s. (B) 
Averaged LFP power spectrum for an external mean rate of 8500 spikes/s 
(solid line) and 5000 spikes/s (dashed line). The black bar delimits the 



gamma-band (30-90 Hz). (C) Raster plot of 200 neurons for the same 250-ms 
interval. (D) Time-resolved firing-rate of the whole population. (E) The 
corresponding power spectrum for the excitatory (black) and inhibitory (red) 
neurons. 
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FIGURE 4 | Modeled network with a GABAergic decay time constant 

Td = 5 ms. (A) Histogram of instantaneous firing-rates of both excitatory 
(black) and inhibitory (red) neurons. The arrows point at the minimum of the 
distribution separating the fast and slow firing mode. (B) Time trace of the 
external excitatory current (brown lines) and the net synaptic current (green 



lines) acting upon the excitatory neurons. The thick lines are averages over 
the entire excitatory population, whereas the thin lines correspond to the 
current impinging on a single neuron. The asterisks label the spiking times of 
this neuron. The horizontal dashed line marks the spiking threshold, i.e., the 
value of / at the bifurcation of Figure 1A. 
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currents acting upon the excitatory population (which by def- 
inition determines the dynamics of the LFP, see section 2.4). 
Figure 4B shows the time traces of the external excitatory input 
(brown lines) and the net synaptic current (green lines). This 
latter current accounts for both the excitatory component (aris- 
ing from the external and the recurrent excitatory spikes) and 
the inhibitory component (arising from the recurrent inhibitory 
spikes). The thick lines are population averages over the excita- 
tory neurons, whereas the thin lines represent the current values 
impinging locally on a given neuron. As shown in this figure, 
the average external synaptic input (thick brown line) lies well 
above the spiking threshold (horizontal dashed line), whereas 
the average net synaptic current (thick green line) is below the 
spiking threshold, due to the inhibitory flow that counteracts 
the strong external excitation. Therefore, the current impinging 
on a neuron [thin green line, corresponding to the 7 syn term of 
Equation (1)] is typically below threshold. Thus neurons spike 
rarely, only when the excitation-inhibition balance is lost during a 
certain time window, in which the external input brings the neu- 
ron above threshold (see asterisks in Figure 4B). From time to 
time these intervals are long enough for several spikes to occur 
in quick succession, giving rise to periods of high-firing-rate (for 
instance around ~1325 ms). On the other hand, if the excitation- 
inhibition balance is briefly lost, an isolated spike is elicited only 
if inhibition is low enough (as we will see in detail in section 3.5), 
giving rise to low-firing-rate events (for instance at ~ 1230 ms). 
Therefore the combined dynamics of excitation and inhibi- 
tion is the basis of the bimodal distribution of instantaneous 
firing-rates. 

According to the previous discussion, the ratio of fast to slow 
firing events is determined by the characteristic timescale of the 
periods in which the excitation-inhibition balance is lost. In order 
to verify this reasoning, we now analyze the behavior of the net- 
work for a longer value of the decay time constant of the 
GABAergic synapses. Increasing to 30 ms leads to a disap- 
pearance of the gamma rhythm (compare the thick green lines 
in Figures 4B, 5B) (Fisahn et al, 1998; Heistek et al., 2010). 



Since the inhibitory currents are slower, the periods in which 
the excitation-inhibition balance is lost are longer, giving rise to 
an increase in the number of fast spiking events (see asterisks 
in Figure 5B). This is quantified in Figure 5A, which shows that 
high firing-rates are much more frequent than low firing ones. 
Thus the inhibition decay time, x L \, determines the principal indi- 
vidual firing-rate mode (mostly single-spike for fast inhibition 
and bursty for slow inhibition). 

3.4. CODING IMPLICATIONS OF THE FIRING-RATE BIMODALITY 

The type I network oscillating in the gamma range, Figure 4A, 
clearly has a sparse activity, in agreement with the sparse coding 
of sensory inputs (Olshausen and Field, 2004; Wolfe et al., 2010). 
We now ask how the bimodality between slow and fast spiking 
regimes, with a band of forbidden firing-rates, affects the cod- 
ing capabilities of the neuronal network. To address this question, 
we establish how the individual firing-rate depends on the exter- 
nal input to the population, rather than on the internal synaptic 
current (as shown in Figure 1 above). In our case, the external 
input is modeled by a set of spike trains perturbing each neuron, 
mimicking either an external sensory input representing a stimu- 
lus or neuronal activity arriving from other areas. As mentioned 
above, the time course of the instantaneous rate is an Ornstein- 
Uhlenbeck process, equal for all the spike trains (but the specific 
realization of the Poisson process is different for each neuron). Its 
mean value is a measure of the input intensity. 

Figure 6 shows the instantaneous individual firing-rate of the 
excitatory neurons composing the network as a function of the 
mean external firing- rate, averaged over the corresponding ISI for 
which the instantaneous firing-rate is calculated, for three differ- 
ent time scales of the GABAergic synaptic dynamics, all of them 
generating LFP oscillations in the gamma range. All three plots 
clearly reveal the bimodal character of the firing-rate distribution 
described above. Two clouds of data points are clearly discernible, 
corresponding to distinct ISIs. The cloud at the top exhibits a 
clear correlation with the external rate, which provides for a stan- 
dard mechanism of rate coding: the higher the external rate is, 
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FIGURE 5 | Modeled network with x d = 30ms. (A) Histogram of 
instantaneous firing-rates of both excitatory (black) and inhibitory 
(red) neurons. (B) Time trace of the external excitatory current 
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(brown lines) and the net synaptic current (green lines) acting upon 
the excitatory neurons. See caption of Figure 4 for more details on 
this plot. 
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FIGURE 6 | Instantaneous firing-rates of the excitatory neurons across 
20 trials as a function of the mean external firing-rate averaged over 
each corresponding ISI. The mean rate of the external spike train is 
8500 spikes/s and is marked with a vertical dashed line. Three different time 
scales of the GABAergic synaptic dynamics, controlled by the decay time 



the faster the resulting firing-rate of the neurons in the network. 
The lower cloud, on the other hand, is associated with longer 
firing periods. Naturally, the instantaneous external rate of that 
low- firing- rate state converges to the mean value of the external 
rate (marked by a vertical dashed line in the three plots), due to 
the large periods over which that mean is calculated. Accordingly, 
for instantaneous external rates clearly above the mean external 
rate the network only responds with high-frequency firing, fol- 
lowing the/-/ curve of the individual neurons (i.e., operating in 
a rate coding mode). For moderate external rates, on the other 
hand, the distribution of firing-rates is bimodal, with low-firing 
and high-firing events coexisting for the same value of the exter- 
nal rate. This allows the network to encode for the phase of the 
global oscillations, since in each one of the two modes the lock- 
ing of the individual spikes to the gamma rhythm is different, as 
described in section 3.5 below. 

The scatter plots shown in Figure 6 also allow us to further 
investigate the origin of the firing-rate gap evident in those plots 
(and in Figure IB above). We can anticipate that the basis of 



constant \a, are considered: (A) xa = 2.5 ms, (B) xa = 5.0 ms (used in the 
rest of the simulations presented in this study), and (C) x d = 10 ms. The local 
density of firing events within each cloud is represented by a normalized 
color scale (with red corresponding to high values and dark blue to low 
values) in order to emphasize that the slow firing mode predominates. 



that gap is the fact that due to the recovery time of the network 
following an increase of inhibition, firing-rates that match the 
duration of network depression are forbidden. To verify that 
expectation, we compare the behavior of the network for dif- 
ferent values of the decay time constant of the GABAergic 
synapses, while keeping the IPSP amplitude and AM PA synapses 
constant. Our results show that longer inhibition [growing from 
panel (A) to panel (C) in Figure 6] increases the range of for- 
bidden firing periods, because the probability of spiking after 
an inhibitory barrage is lower during a longer time interval. 
Accordingly, as increases the low-firing mode reduces its area 
because the slower inhibitory synapses forbid the shortest ISIs 
within this mode. Therefore, the low-firing state is the result of 
a competition between the external slow fluctuations and the 
recurrent inhibition, whose oscillating frequency decreases with 
increasing GABA synaptic decay time. This leads to the counterin- 
tuitive discontinuity in the firing-rate of type I neurons described 
above. In summary, neurons behave like all-or-none detectors of 
rapid stimulus fluctuations, with instantaneous firing-rates faster 
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than the LFP gamma peak encoding for the fast dynamics of the 
stimulus. 

3.5. PHASE LOCKING TO THE GAMMA CYCLE 

We have seen that, even though the firing-rate of the individ- 
ual neurons is far from being tonic, the probability of firing 
across the network varies rhythmically in time with a frequency 
around 45 Hz. This rhythm is generated by recurrent excitatory 
and inhibitory connections, and is revealed in the LFP and firing- 
rate dynamics. Both the synaptic flow (apparent in the LFP) and 
the spiking activity are mutually interacting, given that a decrease 
in synaptic inhibition triggers an increase in the mean firing activ- 
ity of the network. This causal relationship implies that the peaks 
of the population firing activity and the troughs of the LFP are 
displaced in time within the millisecond range, controlled by the 
synaptic delay. 

To further characterize the effect of the global rhythm on the 
firing activity, we have computed STA across trials of both the 
LFP and the inhibitory synaptic current, 7gaba> acting upon each 
neuron. The LFP carries information about the global synaptic 
activity affecting the excitatory neurons of the population (see 
Equation 2), and is therefore a measure of the global activity of 
the network. In contrast, -Tgaba only accounts for the inhibitory 
synaptic current impinging on an neuron from the firing activity 
of its presynaptic inhibitory neurons, and is thus a local measure. 
We have considered the STA in a time window of 70 ms around 
a spike (50 ms pre-spike and 20 ms post-spike). Furthermore, to 
avoid having a previous spike from the same neuron fall within 
the time window being considered, we have only taken into 
account spikes at least 50 ms apart. When computing the STA 
we have classified the spikes, according to the bimodal distribu- 
tion of instantaneous firing-rates shown in Figure 4A, as slow or 
fast firing events, the latter corresponding to network-induced 
bursts. Within the fast firing mode we selected the first spike of 
each burst, since we are interested in the events leading to burst 



initiation. We now discuss separately the mechanisms underly- 
ing the generation of the slow and fast firing modes and their 
relationship with the phase of the global oscillations. 

3.5.1. Fast firing mode 

As shown by the thin lines in Figure 7 A, the high-frequency 
regime appears when the inhibitory synaptic current Tgaba 
impinging on a neuron (i.e., its local inhibition) is close to zero for 
a relatively large time interval right after f = 0. In other words, the 
presynaptic inhibitory neurons must be silent for a long enough 
period of time after the neuron fires, in order for a burst to be ini- 
tiated. In those conditions, it is more probable that the excitatory 
external input brings the neuron above the spiking threshold for 
a sustained amount of time, giving rise to a fast spiking period. 

For the excitatory neurons this occurs preferentially at a high 
level of global inhibition, i.e., at the maxima of the LFP (thin 
line in the top panel of Figure 7B). These two events, a maxi- 
mum of the global inhibition (i.e., of the LFP) and a minimum 
of the local inhibition (i.e., of |IgabaI)> occur simultaneously, 
since it is when global inhibition is strong that some neurons 
can be locally surrounded by strongly inactivated inhibitory neu- 
rons. During this time interval, those neurons fire in bursts due to 
the external excitation, before all the inhibitory neurons become 
excited and prevent the firing of the rest of the population. A 
scheme of the interplay between local inhibition and global acti- 
vation leading to the fast firing mode is shown in Figure 8B. In 
contrast with the excitatory neurons, in inhibitory neurons the 
initiation of the fast firing mode does not depend so strongly 
on the level of global inhibition (thin line in the bottom panel 
of Figure 7B). We hypothesize that these neurons, which have a 
smaller membrane time constant than the excitatory ones (see 
section 2), react rapidly to the external excitatory fluctuations and 
spike regardless of the level of global inhibition, provided their 
presynaptic inhibitory current is zero (thin line in the bottom 
panel of Figure 7A). 




FIGURE 7 | Spike-triggered averages (STA) from (A) the synaptic inhibitory 
current impinging each individual cell and (B) the LFR Spikes from Sis 
larger than 50 ms are considered. The thick lines correspond to those spikes 
prior to a long ISI (slow firing mode), while the thin lines correspond to the short 



ISIs (fast firing mode). These plots are obtained from a single trial, and the 
frontier between long and short ISIs is set at 17.14 ms (58.31 spikes/s) for the 
excitatory neurons and at 11.49 ms (87.01 spikes/s) for the inhibitory neurons. 
The vertical dashed line is the reference time at which a spike is elicited. 
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FIGURE 8 | Scheme of the network spiking activity underlying slow (A) 
and fast (B) firing modes of a given neuron in the network (represented 
by the gray circle). Vertical thick ticks represent spikes. The black and red 
circles represent excitatory and inhibitory neurons, respectively. For 
simplicity, only contacts to the gray neuron are represented (black arrows). In 
the top panel of (A) (slow firing mode), the reference neuron (in gray) is 
surrounded by presynaptic inhibitory neurons that are starting to fire (i.e., the 



level of /gaba. represented by a horizontal red bar, is moderate). The global 
activity of the network (bottom panel, represented by a white horizontal bar) 
is low and the earliest inhibitory spikes prevent the cells from bursting. In the 
top panel of (B) (fast firing mode) the neuron is surrounded by inactive 
inhibitory neurons (i.e., the level of /gaba ' s l° w K whose silence is driven by 
the high activity of the rest of the inhibitory population (reflected in a high 
level of the LFP). 



3.5.2. Slow firing mode 

The low-frequency regime, on the other hand, takes place 
when the local inhibitory current, /gaba> impinging on the 
neuron is momentarily zero, as shown by the thick line in 
Figure 7 A. The rapid increase of inhibition following a spike 



(at t = 0, dashed line) prevents the neuron from spiking again. 
This happens at low levels of global inhibition, i.e., close to 
the minima of the LFP (thick line in Figure 7B). Immediately 
after a spike the firing activity of the population increases, 
and thus the local inhibition grows as well, increasing the 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



March 2013 | Volume 7 | Article 18 | 10 



Sancristobal et al. 



Bimodal firing during gamma oscillations 



distance of individual neurons to threshold (see the scheme of 
Figure 8 A). 

The STA analysis discussed above shows that the simultane- 
ous occurrence of slow instantaneous rates at both the excitatory 
and inhibitory populations (thick lines in Figure 7B) is only pos- 
sible near the minima of the LFP, whereas fast instantaneous rates 
can only be simultaneously present in both populations at the 
maxima of the LFP (thin lines in Figure 7B). Given this inter- 
action between the individual firing modes and the rhythmic 
network dynamics, one can ask whether a partial representation 
of the input is coded in the spiking timing of cells relative to the 
phase of the population oscillation. We now address this ques- 
tion by proceeding to filter the LFP signal around the gamma 
frequency peak (46.88 ± 5 Hz) and assign an instantaneous phase 
to the LFP time series via the analytical signal approach (Hilbert 
phase) (Le Van Quyen et al., 2001). A histogram of the LFP phase 
values at which the neurons spike is shown in Figure 9, with tt 
corresponding to the LFP troughs (vertical dashed line in the 
figure). 

Spike firings are widely distributed over the 2it cycle of the LFP 
phase because, as seen from the raster plot of Figure 3C, the fir- 
ing activity is noisy and the correlation between neurons is weak, 
due to their low coupling and the absence of a common drive. 
We have performed the spike-LFP phase locking analysis for both 
high-firing (thin lines) and low-firing (thick lines) events, taking 
only into account the first spikes of short and long ISIs, respec- 
tively. Figure 9 shows that the slow firing regime of the inhibitory 
neurons (red thick line) is more sensitive to the LFP phase than 
the fast firing regime, resulting in a more pronounced locking 
to the troughs of the LFP, with a delay >2 ms due to the mean 
synaptic delay. The high-firing mode of the excitatory neurons, 
on the other hand, is phase-locked with no delay to the LFP peaks 



around 0 (note that no synaptic delay is considered in the external 
train of spikes, which as mentioned above controls the high-firing 
regime). 

In section 3.4 we have shown that the fast fluctuations of the 
excitatory external rate induce high instantaneous firing-rates, 
and thus this mode encodes the fast dynamics of the input. In 
this section we have additionally shown that the occurrence of 
spikes relative to the LFP phase depends on the firing mode. In 
particular, the high-frequency mode appears with higher proba- 
bility when inhibition is maximal (the LFP phase is 0) and the 
low-frequency mode appears when inhibition is minimal (at the 
rising LFP phases near Tt). This constitutes a mechanism of phase 
coding that complements the rate-coding mechanism depicted 
in Figure 6, according to which fast external fluctuations were 
unambiguously encoded by a high-frequency bursting activity of 
the neurons. 

4. DISCUSSION 

In this paper we have studied the firing characteristics of a neu- 
ral network composed of type I neurons. From a single-neuron 
perspective, both the / -I and PRC have been widely studied and 
related to the kind of bifurcation that takes the cells from rest to 
tonic firing. On the other hand, it is known that cortical neurons, 
wired among them in a complex manner, do not spike period- 
ically but rather sparsely and randomly. The / -I curve of the 
excitatory and inhibitory neural types helps to determine how 
the synaptic connection strength has to be tuned in order to bal- 
ance excitation with inhibition, in agreement with the firing-rates 
of each population. Moreover, according to their phase response, 
type I neurons tend to fire faster within the network than in 
isolation. 

The recurrent connections between the excitatory and 
inhibitory population lead to rhythmicity in the synaptic current 
impinging on a neuron, composed by a succession of a flood of 
excitation and a flood of inhibition. These oscillations are a collec- 
tive phenomenon arising from the interactions among neurons, 
and their frequency can be slowed down with increasing GABA 
synaptic duration. In that way, gamma oscillations in the firing- 
rate activity and LFP signal arise naturally from the recurrent 
connections between excitatory and inhibitory neurons, when 
the synaptic conductances are adjusted to balance excitation with 
inhibition. 

As long as inhibitory neurons fire more intensely than the 
excitatory ones to compensate for their relative small number, 
the network settles in a non-periodic firing state with a signif- 
icant gamma-band component. At frequencies closer or below 
the gamma LFP peak the external drive is balanced by inhibition, 
which leads to a low-firing mode that arises separated from the 
more "natural" fast-firing mode (given by the single-neuron f- 
I curve) by a gap of quasi-forbidden instantaneous firing-rates. 
This leads to a discontinuity in the firing-rates of the neurons, 
similar to what happens in single type II neurons but that does 
not occur in type I neurons in isolation. Thus, this firing-rate 
discontinuity constitutes an emerging dynamical property of the 
network. 

Although the code strategy used by cortical neurons is still 
a matter of debate, bimodal distributions of instantaneous 
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FIGURE 9 | Distribution of LFP phase values at which neurons spike. 

The plot distinguishes between spikes flanking ISIs with fast (thin lines) and 
slow (thick lines) firing modes, for both excitatory (black) and inhibitory (red) 
neurons. Measures were performed over the same data used in Figure 7. 
The dashed line marks the n phase that corresponds to the LFP minimum. 
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firing-rates are found experimentally in both the auditory cor- 
tex (Shih et al, 2011) and the visual cortex (DeBusk et al., 1997; 
Snider et al., 1998). Other modeling works have also shown 
a bimodal distribution of ISIs (Wang, 1998) using intrinsically 
bursting neurons, in contrast with our case in which the bimodal- 
ity appears naturally from the balanced network as in (Roxin 
et al, 2011). Moreover, in Bereshpolova et al. (2011) the recorded 
single unit activity through different awake brain states, shows 
that during alert periods some particular neural types reduce their 
bursting with respect to the non-alert periods. In agreement with 
our results, the alert state corresponds to higher power in the 
gamma range. 

In conclusion, the balance of an excitatory synaptic current 
by a strong inhibitory current yields to a discontinuity in the 
firing-rate of individual neurons forming a neural network. The 
highest instantaneous rates encode fast fluctuations of the exter- 
nal stimulus, while the spiking times of the network occurring 
at moderate fluctuations of the input with respect to the mean 
encode the phase of the LFP. Therefore neural networks can effi- 
ciently implement two coding strategies: (1) a rate code for the 
fast bursting mode, sensitive to rapid changes in the processing 
of stimuli and (2) a phase code for the slower input fluctuations, 
according to which isolated spikes occur at the troughs of the LFP, 



whereas bursts of frequency higher than the gamma peak of the 
global oscillations arise at the peaks of the LFP. This second fea- 
ture might contribute to an internal cortical representation of the 
input. In summary, these results show that cortical population 
activity depends non-trivially on the dynamical properties of the 
underlying neurons, and that global population measures shape 
the firing dynamics of the constituent cells, allowing for multi- 
ple encoding mechanisms to be implemented in networks with 
balanced excitation and inhibition. 
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